stanEx1
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}normal distribution
ex1.stan
normal distribution
data{
int N;
vector[N] y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~normal(m,s);
}mdl=cmdstan_model('./ex1.stan')
y=rnorm(10,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.52 1.06 1.74 1.86 2.55 4.14par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y') 
data=list(N=length(y),y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -13.6463
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -5.59216 7.78142e-05 3.67445e-05 1 1 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.3 seconds.mle## variable estimate
## lp__ -5.59
## m 1.86
## s 1.06mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -6.57 -6.26 1.03 0.76 -8.58 -5.57 1.00 715 752
## m 1.87 1.87 0.43 0.41 1.17 2.56 1.01 697 714
## s 1.30 1.24 0.36 0.28 0.86 1.93 1.00 1295 960mcmc$metadata()$stan_variables## [1] "lp__" "m" "s"mcmc$metadata()$model_params## [1] "lp__" "m" "s"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -6.57 -6.26 1.03 0.76 -8.58 -5.57 1.00 715 752
## m 1.87 1.87 0.43 0.41 1.17 2.56 1.01 697 714
## s 1.30 1.24 0.36 0.28 0.86 1.93 1.00 1295 960
poisson distribution
y=rpois(20,3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 2.00 3.00 3.05 4.00 7.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=length(y),y=y)ex2-1.stan
in case fit poisson dist. to normal dist.
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> m;
real<lower=0> s;
}
model{
y~normal(m,s);
}
generated quantities{
vector[N] y1;
for(i in 1:N)
y1[i]=normal_rng(m,s);
}# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -34.2029
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 -17.643 0.00347767 0.00136638 0.9986 0.9986 24
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -17.64
## m 3.05
## s 1.47
## y1[1] 1.29
## y1[2] 1.93
## y1[3] 3.15
## y1[4] 2.82
## y1[5] 3.54
## y1[6] 3.69
## y1[7] 2.17
## y1[8] 3.85
## y1[9] 2.87
## y1[10] 2.60
## y1[11] 2.47
## y1[12] 0.96
## y1[13] 5.63
## y1[14] 3.03
## y1[15] 1.84
## y1[16] 3.32
## y1[17] 3.20
## y1[18] 0.70
## y1[19] 3.26
## y1[20] 2.78mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.28 -16.90 1.21 0.81 -19.60 -16.19 1.00 852 949
## m 3.05 3.04 0.38 0.36 2.42 3.67 1.00 1747 1065
## s 1.63 1.59 0.31 0.28 1.21 2.16 1.01 825 802
## y1[1] 3.12 3.10 1.74 1.61 0.37 6.04 1.00 1964 1834
## y1[2] 3.06 3.07 1.72 1.63 0.24 5.96 1.00 1754 1759
## y1[3] 3.08 3.11 1.69 1.55 0.34 5.82 1.00 2093 1875
## y1[4] 3.00 3.06 1.70 1.68 0.16 5.70 1.00 1785 1801
## y1[5] 2.98 2.98 1.69 1.61 0.16 5.68 1.00 2021 1870
## y1[6] 3.09 3.08 1.70 1.63 0.32 5.90 1.00 2071 2027
## y1[7] 3.03 2.95 1.71 1.63 0.34 5.92 1.00 1993 2054
## y1[8] 3.01 3.05 1.70 1.59 0.16 5.86 1.00 1532 1593
## y1[9] 2.99 2.96 1.69 1.61 0.25 5.88 1.00 1952 1931
## y1[10] 3.07 3.10 1.74 1.69 0.23 5.95 1.00 1594 1927
## y1[11] 3.03 3.03 1.72 1.68 0.22 5.79 1.00 2015 1995
## y1[12] 3.08 3.08 1.75 1.73 0.20 5.90 1.00 2074 1890
## y1[13] 3.02 3.05 1.65 1.56 0.28 5.70 1.00 2050 1970
## y1[14] 2.98 2.98 1.70 1.69 0.18 5.78 1.00 2034 1895
## y1[15] 3.08 3.08 1.65 1.59 0.43 5.73 1.00 1897 1916
## y1[16] 3.02 2.98 1.73 1.67 0.26 5.80 1.00 1987 1840
## y1[17] 2.99 3.03 1.66 1.56 0.24 5.66 1.00 2021 1849
## y1[18] 3.07 3.07 1.70 1.67 0.30 5.88 1.00 1857 1793
## y1[19] 3.07 3.05 1.71 1.70 0.29 5.86 1.00 2049 2057
## y1[20] 3.09 3.09 1.72 1.56 0.22 5.91 1.00 1776 1725
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
ex2-2.stan
fit to poisson dist.
data{
int N;
array[N] int y;
}
parameters{
real<lower=0> l;
}
model{
y~poisson(l);
}
generated quantities{
array[N] int y1;
for(i in 1:N)
y1[i]=poisson_rng(l);
}mdl=cmdstan_model('./ex2-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -16.1399
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 7.02364 0.000192372 3.67643e-05 1 1 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 7.02
## l 3.05
## y1[1] 4.00
## y1[2] 2.00
## y1[3] 4.00
## y1[4] 3.00
## y1[5] 2.00
## y1[6] 1.00
## y1[7] 3.00
## y1[8] 1.00
## y1[9] 1.00
## y1[10] 4.00
## y1[11] 7.00
## y1[12] 2.00
## y1[13] 1.00
## y1[14] 1.00
## y1[15] 2.00
## y1[16] 4.00
## y1[17] 4.00
## y1[18] 3.00
## y1[19] 4.00
## y1[20] 0.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 7.66 7.93 0.69 0.30 6.29 8.15 1.00 1032 1449
## l 3.10 3.10 0.39 0.38 2.48 3.73 1.00 817 1208
## y1[1] 3.09 3.00 1.81 1.48 1.00 6.00 1.00 1805 1596
## y1[2] 3.11 3.00 1.82 1.48 0.95 6.00 1.00 1666 1870
## y1[3] 3.11 3.00 1.75 1.48 1.00 6.00 1.00 1922 1633
## y1[4] 3.03 3.00 1.80 1.48 0.00 6.00 1.00 1877 1939
## y1[5] 3.10 3.00 1.81 1.48 0.00 6.00 1.00 1998 1622
## y1[6] 3.11 3.00 1.83 1.48 0.00 6.00 1.00 1932 1788
## y1[7] 3.06 3.00 1.79 1.48 1.00 6.00 1.00 1786 1868
## y1[8] 3.01 3.00 1.78 1.48 0.00 6.00 1.00 1835 1744
## y1[9] 3.12 3.00 1.76 1.48 1.00 6.00 1.00 1917 1960
## y1[10] 3.12 3.00 1.77 1.48 0.00 6.00 1.00 1741 1947
## y1[11] 3.04 3.00 1.75 1.48 0.00 6.00 1.00 1973 1988
## y1[12] 3.04 3.00 1.77 1.48 0.00 6.00 1.00 1951 1832
## y1[13] 3.10 3.00 1.78 1.48 1.00 6.00 1.00 1987 2020
## y1[14] 3.01 3.00 1.80 1.48 0.00 6.00 1.00 2104 1909
## y1[15] 3.14 3.00 1.82 1.48 1.00 6.00 1.00 1889 1495
## y1[16] 3.08 3.00 1.85 1.48 1.00 7.00 1.00 1934 1799
## y1[17] 3.10 3.00 1.79 1.48 1.00 6.00 1.00 1928 1994
## y1[18] 3.08 3.00 1.78 1.48 1.00 6.00 1.00 2085 1995
## y1[19] 3.14 3.00 1.79 1.48 1.00 6.00 1.00 1730 1877
## y1[20] 3.12 3.00 1.81 1.48 0.00 6.00 1.00 1656 1811
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
difference test
ex3.stan
mean difference of 2 normal distributions
data{
int N;
vector[N] a;
vector[N] b;
}
parameters{
real ma;
real<lower=0> sa;
real mb;
real<lower=0> sb;
}
model{
a~normal(ma,sa);
b~normal(mb,sb);
}
generated quantities{
real d;
d=mb-ma;
}n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)
mdl=cmdstan_model('./ex3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -29116.3
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 -32.5812 0.000187185 0.000717354 1 1 52
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -32.58
## ma 10.14
## sa 1.08
## mb 10.71
## sb 1.74
## d 0.56mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -34.01 -33.66 1.49 1.25 -36.93 -32.32 1.00 772 1202
## ma 10.13 10.13 0.26 0.26 9.72 10.54 1.00 1290 1090
## sa 1.19 1.16 0.22 0.20 0.89 1.58 1.00 2567 1394
## mb 10.73 10.73 0.43 0.41 10.02 11.44 1.01 824 769
## sb 1.91 1.86 0.34 0.31 1.45 2.56 1.01 1953 1444
## d 0.59 0.59 0.50 0.50 -0.26 1.42 1.01 899 938
d=mcmc$draws('d')
d## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
##
## chain
## iteration 1 2 3 4
## 1 0.18 -0.23 0.96 0.22
## 2 0.30 1.16 0.96 0.58
## 3 0.62 1.09 1.10 0.53
## 4 0.50 0.97 0.91 0.66
## 5 0.52 1.08 0.95 1.13
##
## # ... with 495 more iterationsmean(d>0)## [1] 0.883normal regression
ex4-1.stan
single normal regression
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)
data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex4-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -824269
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpSt0wv7/model-a38248be373.stan', line 14, column 2 to column 22)
## Error evaluating model log probability: Non-finite gradient.
## 42 -53.5663 0.00362798 0.000931153 0.9809 0.9809 108
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -53.57
## b0 -2.64
## b1 3.20
## s 8.83mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -53.06 -52.69 1.42 1.12 -55.83 -51.55 1.01 319 303
## b0 -2.25 -2.61 6.21 5.68 -12.11 8.26 1.02 175 104
## b1 3.20 3.20 0.09 0.08 3.05 3.33 1.02 198 149
## s 10.08 9.81 1.83 1.63 7.55 13.46 1.00 814 917
b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)
ex4-2.stan
single normal regression, prediction from new explanatory variable
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N1] m;
vector[N1] y1;
for(i in 1:N1){
m[i]=b0+b1*x1[i];
y1[i]=normal_rng(m[i],s);
}
}# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex4-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -18297.4
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 33 -53.5663 0.0672733 0.00109384 1 1 75
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -53.57
## b0 -2.64
## b1 3.20
## s 8.83
## m[1] -1.79
## m[2] 33.55
## m[3] 68.89
## m[4] 104.23
## m[5] 139.56
## m[6] 174.90
## m[7] 210.24
## m[8] 245.58
## m[9] 280.92
## m[10] 316.26
## y1[1] -1.50
## y1[2] 33.50
## y1[3] 71.09
## y1[4] 108.73
## y1[5] 160.00
## y1[6] 168.09
## y1[7] 202.29
## y1[8] 248.71
## y1[9] 272.47
## y1[10] 325.25mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.3 seconds.
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -53.04 -52.62 1.38 1.10 -55.70 -51.55 1.01 330 423
## b0 -1.62 -1.96 5.94 5.56 -11.02 9.34 1.02 199 164
## b1 3.19 3.20 0.08 0.08 3.05 3.33 1.02 226 273
## s 9.98 9.73 1.81 1.71 7.53 13.34 1.00 864 787
## m[1] -0.78 -1.12 5.92 5.54 -10.14 10.15 1.02 199 164
## m[2] 34.42 34.14 5.08 4.67 26.41 43.77 1.02 201 125
## m[3] 69.62 69.42 4.28 3.95 62.88 77.37 1.03 206 177
## m[4] 104.82 104.74 3.54 3.31 99.32 111.07 1.02 218 226
## m[5] 140.02 139.98 2.91 2.77 135.50 145.01 1.02 256 363
## m[6] 175.22 175.17 2.46 2.31 171.31 179.37 1.00 396 532
## m[7] 210.43 210.42 2.31 2.20 206.67 214.39 1.01 1106 1062
## m[8] 245.63 245.66 2.51 2.42 241.53 249.81 1.01 2122 1123
## m[9] 280.83 280.87 2.99 2.87 275.75 285.63 1.01 1254 1084
## m[10] 316.03 316.11 3.65 3.63 309.75 321.87 1.01 688 1012
## y1[1] -0.73 -0.64 11.66 11.82 -19.35 19.11 1.00 613 1204
## y1[2] 34.33 33.90 11.24 10.84 15.96 53.16 1.00 670 1026
## y1[3] 69.82 69.86 11.16 10.59 50.96 87.72 1.00 978 1374
## y1[4] 105.03 105.14 10.77 9.97 87.33 122.45 1.00 865 1380
## y1[5] 139.98 140.20 10.75 10.60 122.25 157.46 1.00 1722 1549
## y1[6] 175.51 175.70 10.56 10.32 158.27 192.08 1.00 1797 1835
## y1[7] 210.27 210.08 10.40 10.20 193.35 227.49 1.00 1970 1869
## y1[8] 245.72 245.83 10.12 10.16 229.02 261.62 1.00 1910 1870
## y1[9] 281.00 280.89 10.58 9.76 262.83 298.58 1.00 1798 1646
## y1[10] 315.80 315.58 10.91 10.72 297.93 334.03 1.00 1767 1525
m=mcmc$draws('m')
summary(m)## # A tibble: 10 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 m[1] -0.779 -1.12 5.92 5.54 -10.1 10.2 1.02 200. 164.
## 2 m[2] 34.4 34.1 5.08 4.67 26.4 43.8 1.02 202. 126.
## 3 m[3] 69.6 69.4 4.28 3.95 62.9 77.4 1.03 206. 177.
## 4 m[4] 105. 105. 3.54 3.31 99.3 111. 1.02 218. 226.
## 5 m[5] 140. 140. 2.91 2.77 135. 145. 1.02 257. 364.
## 6 m[6] 175. 175. 2.46 2.31 171. 179. 1.00 396. 533.
## 7 m[7] 210. 210. 2.31 2.20 207. 214. 1.01 1106. 1062.
## 8 m[8] 246. 246. 2.51 2.42 242. 250. 1.01 2122. 1123.
## 9 m[9] 281. 281. 2.99 2.87 276. 286. 1.01 1255. 1085.
## 10 m[10] 316. 316. 3.65 3.63 310. 322. 1.01 688. 1013.smm=summary(m)
y1=mcmc$draws('y1')
summary(y1)## # A tibble: 10 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 y1[1] -0.730 -0.645 11.7 11.8 -19.4 19.1 1.00 614. 1204.
## 2 y1[2] 34.3 33.9 11.2 10.8 16.0 53.2 1.00 670. 1026.
## 3 y1[3] 69.8 69.9 11.2 10.6 51.0 87.7 1.00 979. 1375.
## 4 y1[4] 105. 105. 10.8 9.97 87.3 122. 1.00 865. 1381.
## 5 y1[5] 140. 140. 10.8 10.6 122. 157. 1.00 1722. 1550.
## 6 y1[6] 176. 176. 10.6 10.3 158. 192. 0.999 1797. 1836.
## 7 y1[7] 210. 210. 10.4 10.2 193. 227. 1.00 1971. 1870.
## 8 y1[8] 246. 246. 10.1 10.2 229. 262. 1.00 1911. 1871.
## 9 y1[9] 281. 281. 10.6 9.76 263. 299. 1.00 1798. 1647.
## 10 y1[10] 316. 316. 10.9 10.7 298. 334. 1.00 1768. 1526.smy=summary(y1)
xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)
par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=m),xy,col='black')+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')
ex4-3.stan
single normal regression, prediction from original explanatory variable
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m;
vector[N] y1;
for(i in 1:N){
m[i]=b0+b1*x[i];
y1[i]=normal_rng(m[i],s);
}
}data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex4-3.stan')
mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -53.04 -52.68 1.33 1.10 -55.64 -51.58 1.01 364 350
## b0 -2.95 -2.95 6.14 5.92 -12.39 7.56 1.03 151 186
## b1 3.21 3.21 0.09 0.08 3.07 3.35 1.03 174 315
## s 10.11 9.92 1.76 1.67 7.62 13.49 1.01 781 690
## m[1] 188.92 188.86 2.35 2.26 185.04 192.76 1.01 608 688
## m[2] 309.72 309.62 3.63 3.54 304.00 315.64 1.01 619 1248
## m[3] 175.51 175.48 2.45 2.27 171.50 179.59 1.01 382 650
## m[4] 262.90 262.89 2.76 2.70 258.47 267.47 1.01 1720 1664
## m[5] 176.45 176.41 2.44 2.26 172.46 180.52 1.01 393 649
## m[6] 71.79 71.74 4.32 4.01 65.10 79.27 1.03 153 204
## m[7] 121.32 121.24 3.26 2.98 116.25 126.91 1.03 171 278
## m[8] 110.26 110.20 3.48 3.20 104.84 116.27 1.03 163 251
## m[9] 111.67 111.61 3.45 3.17 106.29 117.62 1.03 163 248
## m[10] 239.56 239.51 2.47 2.49 235.59 243.72 1.01 2227 1609
## m[11] 194.98 194.93 2.32 2.25 191.16 198.77 1.01 722 696
## m[12] 306.21 306.10 3.55 3.46 300.61 312.01 1.01 651 1244
## m[13] 309.30 309.20 3.62 3.53 303.59 315.20 1.01 623 1248
## m[14] 219.51 219.51 2.33 2.34 215.68 223.28 1.00 1694 1571
## m[15] -2.10 -2.10 6.12 5.88 -11.49 8.38 1.03 151 186
## m[16] 153.46 153.45 2.71 2.50 149.05 158.06 1.02 232 426
## m[17] 308.54 308.45 3.60 3.51 302.85 314.43 1.01 630 1248
## m[18] 316.48 316.35 3.77 3.68 310.49 322.70 1.01 569 1102
## m[19] 265.53 265.51 2.80 2.75 261.03 270.14 1.01 1650 1690
## m[20] 179.65 179.59 2.41 2.24 175.71 183.68 1.01 438 665
## y1[1] 188.67 188.62 10.20 9.84 172.14 205.54 1.00 1622 1691
## y1[2] 309.93 310.02 10.94 10.22 290.90 328.43 1.00 1737 1957
## y1[3] 175.71 175.81 10.27 9.70 158.54 192.63 1.00 1869 1870
## y1[4] 262.53 262.67 10.49 10.21 244.88 279.56 1.00 1993 1880
## y1[5] 176.65 176.59 10.80 10.64 159.09 193.76 1.00 1902 1972
## y1[6] 71.86 71.81 11.19 10.92 53.75 90.43 1.00 1008 903
## y1[7] 121.49 121.32 10.83 10.27 103.52 139.28 1.00 1499 1463
## y1[8] 110.72 110.81 10.70 10.34 93.04 128.44 1.00 1164 1599
## y1[9] 111.97 112.07 10.80 10.25 94.21 129.81 1.00 1240 1733
## y1[10] 239.66 239.61 10.69 10.49 222.17 257.29 1.00 2093 2002
## y1[11] 195.09 195.16 10.35 10.06 178.41 212.12 1.00 1949 1723
## y1[12] 306.11 306.36 10.61 10.73 288.36 323.52 1.00 1884 1928
## y1[13] 309.27 309.20 10.94 10.49 291.03 327.03 1.00 1539 1958
## y1[14] 219.16 219.28 10.54 9.79 201.85 236.33 1.00 1791 1968
## y1[15] -2.49 -2.31 12.05 11.36 -22.20 17.87 1.01 647 1304
## y1[16] 153.57 153.71 10.65 10.39 136.56 170.59 1.00 1736 1773
## y1[17] 308.13 308.30 10.78 10.34 290.42 325.66 1.00 1770 1999
## y1[18] 316.78 316.78 10.85 10.75 298.57 333.96 1.00 1432 1874
## y1[19] 265.50 265.31 10.57 10.42 248.04 282.84 1.00 1980 1886
## y1[20] 179.84 179.95 10.58 10.43 162.19 196.68 1.00 1675 1883
y1=mcmc$draws('y1')
smy=summary(y1)
par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
estimate correlation
ex4-41.stan
use covariance matrix and multinormal distribution
data{
int N;
array[N] vector[2] y;
}
parameters{
vector[2] m;
vector<lower=0>[2] s;
real<lower=-1,upper=1> r;
}
transformed parameters{
cov_matrix[2] cv;
cv[1,1]=s[1]^2;
cv[2,2]=s[2]^2;
cv[1,2]=r*s[1]*s[2];
cv[2,1]=r*s[1]*s[2];
}
model{
y~multi_normal(m,cv);
}ex4-42.stan
use covariance matrix and wishart distribution
data{
int N;
matrix[2,2] dp;
}
parameters{
vector<lower=0>[2] s;
real<lower=-1,upper=1> r;
}
transformed parameters{
cov_matrix[2] cv;
cv[1,1]=s[1]^2;
cv[2,2]=s[2]^2;
cv[1,2]=r*s[1]*s[2];
cv[2,1]=r*s[1]*s[2];
}
model{
dp~wishart(N-1,cv);
}ex4-43.stan
use correlation matrix and wishart distribution
data{
int N;
int M;
matrix[M,M] dp;
}
parameters{
vector<lower=0>[M] s;
corr_matrix[M] r;
}
transformed parameters{
cov_matrix[M] cv;
cv=quad_form_diag(r,s);
}
model{
dp~wishart(N-1,cv);
}ex4-44.stan use covariance matrix and multinormal distribution
data{
int N;
int K;
array[N] vector[K] y;
}
parameters{
vector[K] m;
cov_matrix[K] cov;
}
model{
y~multi_normal(m,cov);
}library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)## [,1] [,2]
## [1,] 0.997 0.590
## [2,] 0.590 0.843cor(y)## [,1] [,2]
## [1,] 1.000 0.643
## [2,] 0.643 1.000plot(y)
#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')
data=list(N=n,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -146.414
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 -11.9003 3.45688e-05 3.28354e-05 1 1 22
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -11.90
## m[1] -0.09
## m[2] -0.01
## s[1] 0.97
## s[2] 0.89
## r 0.64
## cv[1,1] 0.95
## cv[2,1] 0.56
## cv[1,2] 0.56
## cv[2,2] 0.80mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.86 -15.50 1.63 1.50 -19.04 -13.84 1.01 539 1060
## m[1] -0.09 -0.09 0.24 0.23 -0.49 0.31 1.00 1413 1004
## m[2] -0.02 -0.02 0.22 0.22 -0.37 0.34 1.00 1450 1105
## s[1] 1.08 1.05 0.19 0.18 0.81 1.44 1.00 1638 1471
## s[2] 0.99 0.97 0.17 0.15 0.74 1.31 1.00 1655 1570
## r 0.59 0.61 0.15 0.14 0.29 0.80 1.00 592 761
## cv[1,1] 1.19 1.11 0.43 0.38 0.66 2.07 1.00 1638 1471
## cv[2,1] 0.65 0.59 0.31 0.26 0.24 1.24 1.00 718 895
## cv[1,2] 0.65 0.59 0.31 0.26 0.24 1.24 1.00 718 895
## cv[2,2] 1.00 0.93 0.36 0.29 0.55 1.71 1.00 1655 1570
#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n
mdl=cmdstan_model('./ex4-42.stan')
data=list(N=n,dp=cov(y)*(n-1))
mle=mdl$optimize(data=data)## Initial log joint probability = -691.138
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 15 -12.2798 0.000110696 0.000833983 1 1 18
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -12.28
## s[1] 1.00
## s[2] 0.92
## r 0.64
## cv[1,1] 1.00
## cv[2,1] 0.59
## cv[1,2] 0.59
## cv[2,2] 0.84mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.05 -14.71 1.22 0.95 -17.51 -13.74 1.00 920 1093
## s[1] 1.08 1.06 0.19 0.17 0.82 1.45 1.00 1351 1384
## s[2] 0.99 0.97 0.17 0.16 0.75 1.30 1.00 1228 1350
## r 0.60 0.62 0.15 0.14 0.34 0.80 1.01 469 713
## cv[1,1] 1.21 1.12 0.46 0.36 0.68 2.09 1.00 1351 1384
## cv[2,1] 0.67 0.61 0.32 0.28 0.28 1.27 1.01 536 930
## cv[1,2] 0.67 0.61 0.32 0.28 0.28 1.27 1.01 536 930
## cv[2,2] 1.01 0.94 0.37 0.30 0.56 1.70 1.00 1229 1350
#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')
m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))
mle=mdl$optimize(data=data)## Initial log joint probability = -50.0993
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 14 -12.2798 0.000249903 0.000125081 1 1 24
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -12.28
## s[1] 1.00
## s[2] 0.92
## r[1,1] 1.00
## r[2,1] 0.64
## r[1,2] 0.64
## r[2,2] 1.00
## cv[1,1] 1.00
## cv[2,1] 0.59
## cv[1,2] 0.59
## cv[2,2] 0.84mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.37 -14.02 1.27 0.97 -16.86 -13.04 1.00 796 786
## s[1] 1.08 1.05 0.18 0.17 0.82 1.41 1.00 1360 1029
## s[2] 0.99 0.97 0.17 0.16 0.76 1.32 1.00 1355 1211
## r[1,1] 1.00 1.00 0.00 0.00 1.00 1.00 NA NA NA
## r[2,1] 0.59 0.61 0.15 0.14 0.30 0.80 1.00 961 1029
## r[1,2] 0.59 0.61 0.15 0.14 0.30 0.80 1.00 961 1029
## r[2,2] 1.00 1.00 0.00 0.00 1.00 1.00 NA NA NA
## cv[1,1] 1.19 1.10 0.43 0.36 0.68 1.98 1.00 1360 1029
## cv[2,1] 0.66 0.60 0.31 0.25 0.26 1.25 1.00 958 816
## cv[1,2] 0.66 0.60 0.31 0.25 0.26 1.25 1.00 958 816
## cv[2,2] 1.02 0.95 0.38 0.32 0.58 1.74 1.00 1355 1211
mdl=cmdstan_model('./ex4-44.stan')
k=2
data=list(N=n,K=k,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -64.0672
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 14 -11.9003 4.65303e-05 0.000682734 1 1 26
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -11.90
## m[1] -0.09
## m[2] -0.01
## cov[1,1] 0.95
## cov[2,1] 0.56
## cov[1,2] 0.56
## cov[2,2] 0.80mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -13.93 -13.52 1.83 1.58 -17.36 -11.75 1.00 561 792
## m[1] -0.09 -0.09 0.28 0.26 -0.54 0.35 1.01 839 806
## m[2] -0.02 -0.03 0.24 0.23 -0.41 0.39 1.01 904 1006
## cov[1,1] 1.46 1.31 0.62 0.49 0.74 2.62 1.00 1262 1064
## cov[2,1] 0.87 0.75 0.49 0.36 0.32 1.81 1.01 1185 916
## cov[1,2] 0.87 0.75 0.49 0.36 0.32 1.81 1.01 1185 916
## cov[2,2] 1.22 1.09 0.53 0.40 0.64 2.21 1.00 1269 990
